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Abstract 

Nonlinear ideal magnetohydrodynamic (MHD) simulations of the propagation and expansion 
of a magnetic "bubble" plasma into a lower density, weakly-magnetized background plasma are 
presented. These simulations mimic the geometry and parameters of the Plasma Bubble Expansion 
Experiment (PBEX) [A. G. Lynn, Y. Zhang, S. C Hsu, H. Li, W. Liu, M. Gilmore, and C 
Watts, Bull. Amer. Phys. Soc. 52, 53 (2007)], which is studying magnetic bubble expansion as a 
model for extra- galactic radio lobes. The simulations predict several key features of the bubble 
evolution. First, the direction of bubble expansion depends on the ratio of the bubble toroidal to 
poloidal magnetic field, with a higher ratio leading to expansion predominantly in the direction 
of propagation and a lower ratio leading to expansion predominantly normal to the direction of 
propagation. Second, an MHD shock and a trailing slow-mode compressible MHD wavefront are 
formed ahead of the bubble as it propagates into the background plasma. Third, the bubble 
expansion and propagation develop asymmetries about its propagation axis due to reconnection 
facilitated by numerical resistivity and to inhomogeneous angular momentum transport mainly 
due to the background magnetic field. These results will help guide the initial experiments and 
diagnostic measurements on PBEX. 
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I. INTRODUCTION 



Cavities with diameters from several to tens of kiloparsecs have been observed in the X-ray 
emission from nearly two dozen galaxies, groups, and clusters.- These cavities are filled with 
magnetic fields and relativistic plasmas that radiate in radio emission from "radio lobes. "-^ 
These observations suggest that the X-ray cavities have formed by shoveling aside thermal 
cluster plasmas by radio-emitting plasmas emanating from galaxies. The interactions of 
radio-emitting outflows with X-ray emitting cluster plasmas lead to shocks, which are a 
candidate for heating cluster plasmas to ^ 1 keV.->2 Past theoretical models of these systems 
assume that such outflows are kinetic energy dominated (so-called kinetic energy dominated 
regime) However, recent observations show that both cluster and radio lobe plasmas 

have appreciable magnetic energy.-^^^ This has led to new models in which radio lobes 
are thought to be gigantic "relaxed" plasmas with kilo-to-megaparsec scale jets providing a 
source of magnetic energy and helicity from the galaxy to the lobes.— ^ However, the details 
of how radio lobe magnetic energy and helicity evolve and interact with the intergalactic 
medium are not well understood.— dkiULH These details depend on underlying nonlinear 
plasma physics, including magnetic relaxation of radio lobe plasmas as they expand against 
a background plasma while being driven by jets, heating of the lobe and background plasmas 
due to reconnection and shocks, and angular momentum transport within the lobe and 
between the lobe and background. 

In order to develop further insights into extragalactic radio lobes, a laboratory plasma 
experiment called the Plasma Bubble Expansion Experiment (PBEX)^ has been built to 
address some of the underlying nonlinear plasma physics issues upon which leading radio 
lobe models are based. The experiment will study the related model problem of a magnetic 
plasma "bubble" relaxing and expanding into a lower pressure weakly-magnetized back- 
ground plasma. A new pulsed coaxial gun will form and inject magnetized plasma bubbles 
(i.e., the lobe) into a background plasma (i.e., the intergalactic medium) formed by a heli- 
con and/or hot cathode source on the HELCAT facility. 19 Experimental parameters can be 
adjusted so that important dimensionless parameters, such as plasma /3, are relevant to the 
astrophysical context. 

Numerical modeling helps guide the experiments and aids the data interpretation. In this 
paper we report initial nonlinear simulation results performed with a new three-dimensional 
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(3D) ideal MHD package,-^ which is a time-explicit, compressible, ideal MHD parallel 3D 
code, using high-order Godunov-type finite- volume numerical methods, in Cartesian coordi- 
nates (x,y,z). The simulations mimic PBEX and use experimentally measured or inferred 
parameters (see Tabled]). 

This paper is organized as follows. In Sec. [TTJ, we outline the problem setup including 
initialization of the bubble and background plasma column. We present the simulation 
results in Sec. 1111] and discussions and implications of our results for the experiment are 
given in Sec. [TV] 



II. PROBLEM SETUP 



In the simulations and experiment, a high density magnetized rotating bubble plasma is 
injected radially into a cylindrical plasma volume with a background magnetic field, as shown 
in Fig. [TJ The injected magnetic configuration is not force-free so that Lorentz forces cause 
the bubble to expand while traveling through and interacting with the background plasma. 
The basic model assumptions and numerical treatments we adopt here are essentially the 
same as those in Li et al.— The nonlinear system of time-dependent ideal MHD equations 
in 3D Cartesian coordinates (x, y, z) is given here: 

§£ + V.(p0) = O, (1) 

^ + V-(pvv + (p + ^)I-BB} =0, (2) 

f>F r / b 2 \ 

— + V- lE + p+ — )v-B(v-B) =0, (3) 



E + P+— Jv-B(v-B) 

— -Vx (vxB) = 0, (4) 

in which p, p, v, B and E are the density, (gas) pressure, flow velocity, magnetic field, 
and total energy, respectively. I is the unit diagonal tensor. The total energy is E = 
p/(7 — 1) + pv 2 /2 + B 2 /2, where 7 = 5/3 is the ratio of the specific heats. Note that a 
factor of V^Lt has been absorbed into the scaling for both the magnetic field B and current 
density j. More details are given in Li et al.— All simulations are performed on the parallel 
Linux clusters at Los Alamos National Laboratory. It should be noted that the details of 
effects such as reconnection and heat evolution cannot be addressed accurately due to the 
ideal MHD model and the use of a simplified energy equation. 



Physical quantities are normalized by the characteristic system length scale Rq, density po, 
and velocity C s o based on the measured or expected values from PBEX. The normalization 
factors are summarized in Table [B Normalized variables are used hereafter. 



A. Background plasma equilibria 

A higher pressure magnetic plasma bubble with spherical radius r& = 1, centered initially 
at Xb — 0, yi, — and Zb = —7.33, is injected along the z axis into a lower pressure 
background plasma with injection velocity v- m j (see Fig. [I]). The stationary background 
plasma is composed of a cylindrical plasma column with radius r p = 6.67 confined by a 
background magnetic field B Xj o(r), where r = \/ y 2 + z 2 . Although PBEX will offer a choice 
of gas combinations for the bubble and background, the initial experiments will likely use 
argon for both, and therefore the initial simulations are based on argon with atomic mass 
of 39.948. 

Force balance of the background plasma along the r direction gives 



For r > r p , B x (t = 0) = B x>0 (r p ), where B Xj0 (r p ) is taken to be 3.65 (75 G), while for 
r < r p , B x (t = 0) is determined by the initial pressure profile: 



The background plasma number density and temperature profiles are given by the fol- 
lowing functions: 



which are a good fit to actual experimental data taken by a Langmuir probe on PBEX, 
where the typical 7 P and jt are taken to be 0.33 and 0.14, respectively. Thus, the initial 
pressure profile of the background plasma is p(r, t = 0) = n p T p \ t=0 oc exp[(7 T — j p )r]. 

B. Magnetized bubble plasma 



p(r) + 



B x , (r) 2 _ B xfi (r p ) 2 



2 2 




n p (r,t = 0) = 1.06 exp(— 7 p r), T p (r,t = 0) = 2.7exp(7 T r) . 



A higher pressure magnetized plasma bubble is generated and injected by a coax- 
ial gun source. It is well established empirically in coaxial gun spheromak experiments 
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that, under proper conditions, a spheromak "magnetic bubble" will be formed by the gun 
discharge.— In the simulations reported here, the bubble structure is similar to the 
one given in Li et al.— 

The number density profile of the bubble plasma with radius = 1 is given by 



n b oc r 2 c exp[-r 2 - (z c - z b ) 2 } 



2n 

) 



up to a normalization coefficient n bQ = 100 and a uniform temperature T b0 = 10, where 
r c = \J x 2 + y 2 and z c = z (see Fig. CQ. The density profile used here has its peak shifted 
from the center of the bubble, approximating a spheromak, and is therefore different from 
the uniform density profile used in Li et al.— 

The bubble magnetic field is determined by three key quantities: the length scale of the 
bubble magnetic field r# = 1, the amount of poloidal flux and the index a, which is the 
ratio of the bubble toroidal to poloidal magnetic fields. For simplicity, the bubble magnetic 
field -Bbubbie is also assumed to be axisymmetric. The poloidal flux function ty p is specified 
as: 

^ p cc r 2 c exp[-r c 2 - (z c - z b ) 2 } . (5) 
The poloidal fields, up to a normalization coefficient Bbo = 48.7 (1000 G), are: 

-£>bubble,r c — ~ , £>bubble,2 c — ^ , ID J 

r c oz c r c or c 

while the toroidal magnetic field is 

^bubbie.^c = = ar c exp[-r 2 - (z c - z b ) 2 } . (7) 

' c 

The azimuthal component of the bubble Lorentz force is zero, but the total azimuthal Lorentz 
force due to the combined fields and currents of the bubble and the background plasma may 
be non-zero. 

The bubble also has uniform injection velocity v- m j and uniform rotation angular speed 
Vt = V^7tVA,o/ r bi where Va$ = B b0 / yJAnp b0 = 4.87. Please note that this is a strong 
rotation, possibly having strong influence on the stability of the bubble (Sec JIII B~2j) and the 
expansion of the bubble in the x-y plane (Sec JIII BT) . 
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C. Computational domain 

The total computational domain is \x\ < 9, \y\ < 9, and \z\ < 9, corresponding to a 
(54 cm) 3 box in actual length units. The numerical resolution used here is 400 x 400 x 400, 
where the grid points are assigned uniformly in the x—, y—, and z— directions. A cell Sx 
(= 5y = 5z = 0.045) corresponds to 0.135 cm. We use "outflow" boundary conditions at 
every boundary, i.e., setting all values of variables in the ghost zones equal to the values in 
the corresponding active zones, which is the simplest approach possible. This technique is 
accurate for supersonic outflow but not for subsonic outflow. This simplified boundary con- 
dition limits our ability to predict the transit time, the time for the bubble to travel through 
the background plasma, and to study the detachment problem, i.e., under what conditions 
the bubble would separate from the wall boundary. More accurate boundary conditions will 
be implemented in future work. Here, we focus on the interaction of the bubble plasma with 
the background plasma before the structures have reached the boundaries. 

III. SIMULATION RESULTS 

In this section we present ideal MHD simulation results on the nonlinear evolution of 
a magnetic "bubble" plasma propagating and expanding into a lower pressure background 
plasma. The results are organized into three primary topics: (1) global evolution of the 
bubble-background system and interface, (2) internal bubble evolution, and (3) angular mo- 
mentum transport both outside and inside the "bubble." Key findings include the formation 
of both an MHD shock and a reverse MHD slow-mode wavefront as a result of the bubble 
propagating into the background plasma, and the outward transfer of azimuthal angular 
momentum inside the bubble due to advection and inhomogeneous transport outside the 
bubble due to the background magnetic field. Please note that all physical quantities, such 
as the magnetic field B and flow velocity v, presented in this section are the total value due 
to both the bubble and background plasmas. 

A. Global evolution of the bubble-background system and interface 

In this subsection, we examine the evolution of the global bubble-background system and 
the interface between the two plasmas. 
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1. Bubble propagation and expansion 

Here we discuss the time evolution of the magnetic bubble, showing selected physical 
quantities using 2-D x—z slices at y — 0. The density distributions at various times (t = 
0, 0.25, 0.5, 1.0) are shown in Fig. [2] with a = vTO and injection velocity = 0.18Va,o- I n 
this example a = a/10 corresponds to the bubble having a minimum initial Lorentz force 
(see discussion below). At t = 0.5, we see that the initial peak-shifted high density magnetic 
bubble has been transformed into a "crab" (due to the low a, see discussion below), bounded 
by one MHD shock (see Sec. IIII A 2j) and one reverse slow-mode compressible MHD wavefront 
(Sec. IIII A 3[) . Low-density cavities (a factor of dozens of times of magnitudes smaller than 
the peak density) exist both between the shock and wavefront and in the post- wavefront 
region. At t = 1.5, the shock has reached the other side of the computation domain, while 
the wavefront is located at z ~ 1. The bubble is still in the middle of the background plasma. 
The simulation after t = 1.5 is not accurate due to the simplified boundary conditions. 

The value of a determines the strength of the initial Lorentz force in the bubble and 
consequently how the bubble expands and evolves. We first test the influence of a on the 
bubble evolution with fixed injection speed f; n j. The simulations show that the results are 
insensitive to a except for large a = 15 when the bubble expands more in the direction 
of the injection, leading to a growing "mushroom" (Fig. [3tleft)). With smaller a — 1, the 
bubble expands more transversely to the direction of injection, resulting in a growing "crab" 
(Fig. fright)). This can be understood from the initial poloidal Lorentz force due to the 
coupling of the bubble's current and its own field (Eq. 22 and Eq. 23 of Li et al.— ): larger 
a would have a positive axial Lorentz force but negative radial Lorentz force, resulting in 
collimation, while smaller a would have a negative axial but positive radial Lorentz force, 
resulting in radial expansion. In the experiment, it is expected that the ejected bubble 
will quickly reach a nearly force-free state. Therefore, hereafter, we assume a = vTO, 
corresponding to minimum initial bubble Lorentz force. 

Given reasonably low injection velocity (v- m j ^ Va,o), there are always an MHD shock 
and a reverse slow-mode compressible MHD wavefront, whose structures, evolution, and 
propagation characteristics are essentially the same. And the injection velocity has little 
influence on the shock speed and wavefront speed (see Table HT|) . which implies that the 
shock and wavefront result from the expansion of the bubble due to the Lorentz force, 
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rather than a "piston effect" of the bubble propagation. However, the "piston" effect could 
become important with larger injection speeds (v^ ^ Va,o) (see Table |H| . which, however, 
will not occur in the experiment. For the following analysis, we will focus on the case 
^inj = 0.18Va,o = 0.88 as being representative for the experiment. 

2. Identification of an MHD shock 

The expansion of the magnetic bubble generates a leading MHD shock and a trailing 
reverse slow-mode compressible MHD wavefront. The details of these two structures are 
presented in this and the next subsections. The characteristics of the shock propagation 
are similar in the x— and z— directions. Thus, we will present results and analysis in the 
z— direction only. 

Figure IU with v in j = 0.18Va )0 , displays several physical quantities along a line with 
(x,y) = (0,0) in the ^-direction at t = 0.5. Hereafter we define axial direction as the 
direction along z-axis, x-y plane as toroidal plane and x-z plane as poloidal plane. Several 
features can be identified. First, an MHD shock can be seen around z = —1.035 in the 
profiles of p, B x [see Fig. HJleft)] and V z [see Fig. Upright)]. In the vicinity of x ~ 0, since 
the magnetic field lies in the shock plane and is perpendicular to the shock normal, this 
shock is identified as a perpendicular shock in this region. This MHD shock is a fast shock 
whose properties are very close to an ordinary field-free shock. As shown in Fig. H](left), the 
magnetic field components B y and B z change very little across the shock. 

It is important to verify that this is indeed an MHD shock by comparing the simulation 
results to the the shock jump conditions. Choosing the velocity frame so that the shock 
is at rest (shock velocity is Vg) and simplifying the notation for the problem, we represent 
quantities in the upstream region by a superscript and those in the downstream region 
by no superscript. Then there are 8 known quantities B®, By = 0, B° z = 0, p°, p°, V® = 0, 
Vy = and V® = —Vs. There are 8 unknown quantities in the downstream region, p, p, V x , 
V y , V z , B x , By and B z . Thus we need 8 conditions to specify them. We consider the 1-D 
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ideal MHD shock jump conditions for simplicity. The equations are as follows: 




(10) 
(11) 



(8) 



(9) 



pV v V x -B y B z = 0, 



pv z \(v x 2 + v y 2 + v?) + 




V X B X B Z - V y B y B z + B 2 X V Z + ByV z 




(12) 



B z = B 



■o 



(13) 
(14) 



V z B y -B z V y = 0, 
V Z B X -B Z V X = V Z °B° X . 



(15) 



A MATLAB code was used to solve this nonlinear system of equations, given the values in 
the upstream region: p°, p°, V® and B x . The results of V x , V z , B x and B y matches pretty 
well (see Table IIH]) . The nonzero simulation values of V y and B y result from 3-D effects. 
The relatively large differences seen in the values of p and p are possibly due to nonzero 
numerical diffusion in the simulations. 

3. Reverse slow-mode compressible MHD wavefront 

There is an MHD wavefront at z = —2.835, as seen in both panels of Fig. HJ where B x 
and B y have their local minimum, and p, C s and v x have their local maxima. The nature of 
this MHD wavefront can be identified by plotting the axial pressure profiles along the line 
(x,y) = (0, 0) at t = 0.5, as shown in Fig. [51 A transition occurs near z = —2.835, where 
an increase in gas pressure p is accompanied by a decrease in magnetic pressure p m = B 2 /2. 
The transition is identified as a reverse slow- mode compressible MHD wavefront. Between 
the shock (z = —1.035) and the wavefront (z = —2.835), the magnetic field lines are 
compressed and some thermal energy has been converted into magnetic energy. Therefore the 
gas pressure has an abrupt decrease while the magnetic pressure increases rapidly between 
the shock and wavefront (see Fig. [5]). 
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From Fig. [6j it can be seen that anti-parallel reconnection of the B x component occurs 
around (x,z) = (—5, —4), facilitated by numerical diffusion although the simulation is per- 
formed with an ideal MHD code. It is possible that B y component reconnection takes place 
around (x,z) = (—4,-4), where the bubble field is not exactly anti-parallel to the back- 
ground field (more details in Sec. IIII A 4|) . The reconnection is driven by magnetic field line 
compression due to the reverse slow- mode compressible MHD wavefront. From Fig. [71 two 
strong toroidal current sheets are observed at the reconnection layer and the MHD shock. 
These are also due to the compression of the magnetic field lines due to the wavefront and 
shock, respectively. The reconnection and the shock/ wavefront convert normal velocity into 
tangent velocity and convert kinetic energy into thermal energy. Because this is an ideal 
MHD simulation, the details of the reconnection are not expected to be accurate. We are 
only interested here in the qualitative effects of the bubble evolution due to reconnection. 

4- Force in the z- direction 

The evolution and propagation of the magnetic bubble can be further understood by 
examining the various forces along (x, y) = (0, 0) at t — 0.5, which are displayed in Fig. [HJ 
The MHD shock breaks the initial background equilibrium. The passage of the shock wave 
heats the gas and alters its pressure gradient. The axial flow is pushed forward by both the 
gas pressure gradient and Lorentz force at the MHD shock while it is dragged back behind 
the shock, resulting in an axial deceleration of the gas in the postshock region. Although 
the Lorentz force tries to accelerate the axial flow at the MHD wavefront, the gas pressure 
holds it back. Therefore the MHD shock will be driven forward and eventually separated 
from the wavefront, which leads to a cavity of depleted density between the shock and the 
wavefront. 

Figure M displays axial profiles of various forces along (x,y) = (—4,0) and (x,y) = (4,0) 
at t — 0.5. The locations are chosen to be where anti-parallel reconnection (left panel of 
Fig. [9]) occurs and its reflection about z = (right panel of Fig. [9]). This figure clearly shows 
the difference between the two locations: the Lorentz force changes sign (similar to Fig. [HJ 
on the left hand side while keeping the same sign (negative) on the right hand side, which 
is consistent with reconnection happening on the left hand side and not on the right hand 
side. At both x locations, the toroidal current densities j y are much larger than the poloidal 
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current densities; however, the toroidal field component B y is close to zero. Therefore \j y B x \ 
is much bigger than \j x B y \. Anti-parallel field lines on the left hand side result in a sign- 
change of the axial Lorentz force, while the Lorentz force on the right hand side does not 
change sign since B x does not change sign there. The sign change of the Lorentz force is 
necessary for reconnection since this is the driving force to pull the field lines from either 
side of the current sheet together to reconnect. Also, the total force is more negative on 
the right hand side compared to the left hand side. This means that the axial flow in the 
right hand side is slowed down more quickly than the left hand side, which leads to the 
asymmetry of the shock propagation across the 

The magnetic bubble evolves into a nearly quasi-force-free state (see Fig. [TP]) , with a 
Lorentz force that scales in time roughly as: 

-^Lorentz (t) r t \ 

- exp( 



-^Lorentz 1 1=0 ^relaxation 

The time scale of this relaxation r re iaxation is dependent on the value of a, which determines 
the amplitude of the initial Lorentz force. Larger initial Lorentz force leads to quicker re- 
laxation. For a = 1, r re i axa tion = 0.379; for a = y/lO, r relaxa tion = 0.386; and for a = 15, 
Relaxation = 0.120. The relaxation happens on the order of the Sweet-Parker reconnection 
time rsp = y/r a T res (see definitions and estimates of r a and r res in Sec. HVj) . For example, 
rgp « 0.3 for a = vTO, which is similar to the Lorentz force relaxation time of 0.386. How- 
ever, it should be noted that the Lorentz force in the shock/ wavefront is always significant 

(Fig. ED. 



B. Internal bubble evolution 



In this subsection, we examine the evolution and properties of the bubble itself, including 
a simple kink stability analysis. 



1. Bubble density, velocity and magnetic field evolution 

Density (Fig. [TTT) and fluid velocity vector (Fig. [T2l plots in the x — y plane at different 
times both demonstrate that the initial fast-rotating spheromak-like magnetic bubble evolves 
into a much larger slow-rotating, fast-expanding elliptical structure with maximum density 
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reduced by 20 times, while the density at the wavefront and the shock increases by 8 and 
2 times, respectively. The center of the bubble shifts away from the original propagation 
axis with (x,y) = (0,0). Figure [T51 displays the vector magnetic field plots in the x — y 
plane at different times. From the figure, we can see that the bubble field still keeps a 
spheromak-like configuration and the expansion of the bubble pushes away the background 
fields, thus compressing them. Some reconnection happens in the region y > because the 
bubble toroidal field component is opposite to the background field there. These figures 
show that the background field breaks the symmetry of the system, which is consistent with 
the results of Sec. |niA4|and Sec. IfflCl 

It is worth noting that the initial expansion of the bubble in x-y plane results from 
the non-free initial Lorentz force as well as the centrifugal force due to the strong initial 
rotation of the bubble, although the latter quickly slows down to a small value because 
of the conservation of angular momentum associated with the initial quick expansion and 
the possible Kelvin-Holmhotz instability associated with the initial strong toroidal velocity 
shear. 



2. Bubble stability 

Spheromak-like bubble plasmas are subject to current- driven kink instabilities. Figured 
shows a snapshot of the axial current density j z at t = 0.5. The axial current flow follows a 
semi-closed (it will close outside the out-flowing boundary) circulating path, flowing along 
the central axis (the "forward" current) and returning along the bell-shaped path on the out- 
side (the "return" current).— Fig. [15] shows a snapshot of the configuration of the magnetic 
field B, which indicates that a tightly wound central helix is overlapped with the "for- 
ward" current, and a loosely wound helix is overlapped with the "return" current. Given 
a helical magnetic field, this axial current- carrying cylindrical plasma column is subject 
to a current-driven instability (CDI).— >2£ However, we do not see any visible evidence of 
any current- driven instability in this case (a = vTO). The well-known Kruskal-Shafranov 
criterion2^i2£ for MHD kink instability in cylindrical geometry can be written as: 24 

ff ( fl ) = J^L- < 1 . (16) 

^ J z, total 
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where q is the safety factor, tp p rs ira 2 B z (a) is the total poloidal magnetic flux, i^totai is 
the total axial current, a and L are the column radius and length, and B z is the axial field 
component. Safety factors less than 1 are unstable to the CDI kink mode. The safety factor 
in this case (Fig. [Inland Fig. |T5|) is q(a) ~ 2.8 at t — 0.5, which is bigger than 1. Therefore it 
is expected to be CDI stable, which is consistent with the simulation results. The simulation 
with a = 15 gives q(a) ~ 1 at t = 0.125, which is marginally unstable to CDI according to 
Eq. [16j However we do not find evidence of unstable CDI modes in this case either. "Line 
tying" (important in the experiment, not present in the simulations reported here due to the 
"outflowing" boundary conditions used in this paper )2L2£i22 an( } other stabilization effects 
such as "dynamic relaxation",— internal strong rotation^^ and external gas pressure, etc., 
could raise the stability threshold. A more detailed stability analysis of the magnetic bubble 
is beyond the scope of this paper. 

C. Angular momentum transport 

Since the bubble is rotating about the z— axis initially, the bubble has initial net angular 
momentum. Conservation of azimuthal angular momentum will slow the bubble's rotation 
since some angular momentum will be transported to the background plasma.— It is a key 
nonlinear plasma physics question to address how this angular momentum evolves. 

For an ideal MHD flow, the azimuthal angular momentum conservation equation in cylin- 
drical coordinates (r, ip, z) is:— 



where e v is the unit vector in the azimuthal direction, the p subscript refers to a poloidal 
magnetic-field component (i.e., the r or z component), and B"i = B 2 + B\. There are no 
source terms in this equation, i.e., angular momentum may be redistributed in the fluid but 
never destroyed. The numerical diffusion present in the simulations would transport some 
angular momentum as well. However, the influence of this transport would be highly limited 
in the shock/ wavefront regions and negligible elsewhere. The first term in the bracket Tcpv^v, 
the so-called "advection angular momentum flux" r a d ve ction, is the angular momentum flux 
vector due to the advection, which is defined, in Cartesian coordinates, as: 



0_ 

dt 




3 



0. 



(17) 



r a dvection = p(xv y - yv x )v = p(xv y - yv x )(v x x + v y y + v z z) , 



(18) 
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where x, y and z are the unit vectors in the x—, y— and z— directions respectively. The sec- 
ond term in the bracket —r c B v B p , the so-called "Maxwell angular momentum flux" TMaxweii, 
is the angular momentum flux vector due to the Lorentz force, which is defined in Cartesian 
coordinates as: 

TMaxeii = -{xB v - yB x )B p = -{xB y - yB x )(B r e r + B z e z ) 
= -{xB y - yB x ) [(B x cos 2 9 

+ B v sin 6 cos 9)x + (By sin 2 + B x cos 9 sin 9)y + B z z] , (19) 

where 9 is the polar angle with tan 9 = x/y and e r , e z are the radial and axial unit vectors 
in cylindrical coordinates, respectively. They both contribute to the angular momentum 
transport in every direction. The third term, the so-called "pressure angular momentum 
flux" r pressure , is the angular momentum flux vector due to the effective pressure, which is 
defined in Cartesian coordinates as: 

f pressure = (xf) ~ yx)[ P + (B? + B\) /2] , 

where B r = B x cos9 + B y sin9. This term does not have a z— component, i.e., it only 
distributes azimuthal angular momentum in the toroidal plane (x-y plane). The total angular 
momentum flux r to tai is defined as: 

■Ttotal radvection -^Maxwell ~\~ Tpressurc • 

1. Angular momentum transport in the x-y plane 

Figure QH displays vector plots of r advect ion, r MaX weii, r prcssure and r to tai in the toroidal 
plane, i.e., x — y plane at z — —6, which coincides with the bubble. Pressure simply trans- 
ports angular momentum in an anti-clockwise direction in this plane. Inside the bubble, 
advection transports angular momentum outward, which is due to the expansion of the 
bubble. The Maxwell angular momentum flux only has a radial component in the plane 
(Eq. fT9i) . This flux is dominant outside the bubble since the radial bubble field component 
decreases so quickly that it is much smaller than the background field at t — 0.5 [see also 
Fig. [T3l(right)] due to the quick relaxation (Sec. IIII A 4[) . Interestingly the Lorentz force 
transports angular momentum inward at the top left and bottom right regions and outward 
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at the top right and bottom left regions. Therefore the total effect is to transport angu- 
lar momentum: (1) angular momentum outward inside the bubble; (2) along the negative 
y— axis for x < and along the positive y— axis for x > outside the bubble. The maximum 
angular momentum transport happens at the edge of the bubble. Inside the bubble some 
angular momentum has been transported from the center of the bubble to the edge of the 
bubble, while outside some angular momentum has been transported from the top left to 
the bottom left regions, and from the bottom right to the top right regions [Fig. ITTtnght)]. 
Thus, the uniformly rotating bubble expands and its inner region ceases to rotate and then 
rotates oppositely in the long run, while the neighboring plasma starts to rotate differen- 
tially. The top right and bottom left regions rotate in the same direction as the original 
bubble, while the top left and bottom right regions rotate in the opposite direction (Fig. fTTI) . 
which results in shears. This explains why, between the shock and wavefront, the advection 
transports the angular momentum in negatively when x < while positively when x > 
since the shock is always propagating outward (Eq. [TBI) . 

2. Angular momentum transport in the x-z plane 

Angular momentum transport in the poloidal plane, i.e., the x-z plane at y = 0, is 
presented in Fig. [JH] (r pressure is zero on this plane). Inside the bubble, r a d ve ction due to the 
expansion of the bubble transports angular momentum outward normal to the wavefront, 
while r Maxwe n due to the bubble field redistributes angular momentum inside, transporting 
angular momentum clockwise on the left hand side and anti-closewise on the right hand side. 
This can be understood from the spheromak-like magnetic field configuration of the bubble. 
The total effect is to transport net angular momentum from the right hand side to the left 
hand side of the bubble edge, leading to positive angular momentum on the left hand side 
and negative angular momentum on the right hand side of the bubble edge (Fig. [T9l) . 

IV. SUMMARY & DISCUSSIONS 

In this paper we presented initial nonlinear ideal MHD simulation results of the expansion 
of a magnetic bubble into a lower pressure weakly magnetized background plasma. The sim- 
ulations mimic the ongoing experiment PBEX, except that we use simplified "out-flowing" 
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boundary conditions and ignore collisional effects. A high-density magnetized bubble is 
injected into a cylindrical background plasma. The bubble evolution is dependent on a, 
with larger a resulting in an axially expanding bubble like a growing "mushroom" and 
smaller a producing a "crab-like" shape expanding normal to the direction of propagation. 
The expansion of the bubble generates one leading MHD shock and one trailing reverse 
slow-mode compressible MHD wavefront. The shock/ wavefront speed is independent of the 
injection velocity for injection velocities similar or below Va- In the x-z plane anti-parallel 
reconnection takes place on the left hand side of the wavefront, where the bubble field is 
opposite to the background field, while in the x-y plane the reconnection takes place if y > 0. 
The azimuthal angular momentum is transported outward from the center to the edge of 
the bubble by advection, and from the right hand side to the left hand side of the bubble 
edge by the Maxwell torque. Outside the bubble, angular momentum is transported from 
the top left to the bottom left and from the bottom right to top right by the combination 
of Maxwell and pressure angular momentum fluxes. The initial uniformly rotating bubble 
quickly evolves into a quasi-force-free state, and the center of the bubble ceases to rotate 
and starts to rotate oppositely in the long run, while the outside neighboring plasmas start 
to rotate differentially: the top right and bottom left regions possessing the same rotating 
direction as the original bubbe while the top left and bottom right regions possessing the 
opposite rotating direction, which forms shears in the system. 

The background magnetic field breaks the symmetry of the system along the propagation 
axis. For comparison, simulations with 2 orders of magnitude lower background field (and 
correspondingly lower density and temperature in order to preserve the equilibrium of the 
background plasmas) were also performed. The results show much better symmetry along 
the propagation axis. More detailed studies of the role of the background field will be the 
subject of future work. 

From Table IIVI the resistive dissipation time due to numerical diffusion is inferred to be 
r res ~ 1.2, which is longer than the time (t = 0.5) at which the shock has reached the x 
and y boundaries. The Alfven time can be calculated as r a = L res /VA, r es, where L rcs (~ 0.5) 
is the typical length of the reconnection layer and VA, re s is the Alfven speed (~ 8) at the 
reconnection layer. This gives r a ~ 0.06. Thus the effective Lundquist number ^effective in 
the simulations is around ^effective = T res /r a ~ 20. Please note that this is the Lundquist 
number associated with the reconnection layer and the numerical diffusion used here is 
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the upper limit of the numerical diffusion in the simulations. The estimate of the mean 
experimental Lundquist number ^experiment = (V^Lo/^piasma) is around 200, where L ~ 18 
is the characteristic length of the experimental facility, V® is the initial Alfven speed, ^ p i aS ma 
is the magnetic resistivity of the plasma at the initial state based on Braginskii's formula 
and (...) indicates the volume average. Although there is one order of magnitude difference 
between them, we expect our ideal simulations with numerical diffusion to give a reasonable 
estimate of the physical quantities in the real experiment. 

Another issue important in the experiment is the transit time. Although the value of the 
transit time is somewhat related to the boundary conditions, our simulations with simplified 
boundary condition show that the background plasma column radius (r p = 20 cm) is large 
enough to allow the bubble to relax substantially. From Fig. [2J it is seen that the bubble 
is still inside the background plasma at t — 1.5. However, better boundary conditions are 
needed for the simulations to be meaningful after the shock has reach the boundaries. 

The appearance of the MHD shock and wavefront suggests that our experimental facility 
may provide a unique opportunity to study MHD shocks in a laboratory plasma. However, 
we emphasize that these conclusions are based on ideal simulations (with numerical diffusion) 
and that the boundary conditions are not realistic. This paper is intended as a preliminary 
exploration of PBEX. We have not attempted to model many of the complexities of a realistic 
experiment. In future papers, we will study collisional effects and boundary conditions closer 
to those of the planned experiment; work in progress indicates that these will modify the 
results. 
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Figure 1: Schematic of the simulation geometry of PBEX showing also the coordinate system. In 
the texts, the direction along z— axis is defined as axial direction, x — y plane is defined as the 
toroidal plane while x — z plane is defined as the poloidal plane. 



20 




21 




Figure 3: (color) Density (natural logarithmic scale) for a = 15 (left) and a = 1 (right). 




Figure 4: Axial profiles of physical quantities at (x, y) = (0, 0) and t = 0.5. Left: density p 
and magnetic field components (B x , B y , B z ). Right: sound speed C s and velocity components 

(V X ,Vy,V z ). 
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Figure 5: Axial pressure profiles at (x, y) = (0, 0) and t = 0.5. 
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Fi gure 7: Axial profiles of current density components (jxi jy? jz 

) at t = 0.5. Left: (x,y) = (-4,0). 

Right: (x,y) = (4,0). 
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Figure 8: Axial profiles of various forces at (x, y) = (0, 0) and t = 0.5. F p , the pressure gradient; 
Fj xB , Lorentz force; F tol = F p + F JxB . 
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Figure 9: Axial profiles of various forces at t = 0.5. F p , the pressure gradient; -FjxB, Lorentz force; 
*toi = Fp + i ? JxB- Left: (x,y) = (-4,0). Right: (x,y) = (4,0). 



25 



.000 



s 0.100 



0.010 



0.00 




Figure 10: Maximum absolute value of the axial Lorentz force over the initial value inside the 
bubble versus time. (Dash line: a = 1; dash dot line: a = \/l0; long dash line: a = 15.) 
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Figure 12: (color) Vector plot of flow velocity v at t = 0.5 in the x-y plane. Arrows: flow velocity 
components v x and v y \ color: flow velocity v z . Left: t = 0, z = —7.5. Right: t = 0.5, z = —6. 
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Figure 13: (color) Vector plot of magnetic field B at t = 0.5 in the x-y plane. Arrows: magnetic 
fields B x and B y ; color: magnetic field B z . Left: t = 0, z = —7.5. Right: t = 0.5, z = —6. 
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Figure 14: (color) Axial current density j z at t = 0.5 in the x-z plane at y = 0. 
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Figure 15: (color) Vector plot of magnetic fields B at t = 0.5 in the x-z plane at y = 0. Arrows: 
poloidal magnetic fields B x and B z . Color: toroidal magnetic field B y . 



Physical Quantities 


Description 


Normalized Units 


Typical Values 


R 


Length 


Ro 


3 cm 


V 


Velocity Field 


CgO 


2.0 x 10 5 cm s" 1 


t 


Time 


Ro/Cso 


1.50 x 10~ 5 s 


n 


Number Density 


n 


10 12 cm" 3 


P 


Density 


po = n M a 


6.67 x 10~ n g cm" 3 


P 


Pressure 


PoC% 


2.67 dyn cm -2 


B 


Magnetic Field 




5.79 Gauss 


T 


Temperature 


T 


1 ev 


^Lorentz 


Lorentz Force 


F 


11.2 dyne 



°Mo is the atomic mass. In the experiment it is the argon atomic mass Mq = 6.67 x 10 23 g. 

Table I: Physical quantities, normalization constants, and values. 
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Fi gure 16: Vector plot of the angular momentum fluxes due to advection radvection 

(top left), 

Lorentz force Maxell (top right), pressure r pressure (bottom left) and r Total (bottom right) in the 
x-y plane at z = —6 and t = 0.5. 
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Vinj/VAfl 


Slow-Mode Wave Front Speed 


Perpendicular Shock Speed 


0.0035 


3.6 


7.2 


0.18 


3.6 


7.2 


0.5 


3.6 


7.2 


1.0 


3.6 


7.2 


1.8 


3.6 


10.8 


3.5 


10.8 


18 


35.4 


160 


180 



Table II: Shock speed and wavefront speed vs. injection velocity with a = \/l0. All velocities are 
normalized to the sound speed C s q. Va,,o = 4.87. 
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Fi gure 18: Vector plot of the angular momentum fluxes due to advection radvection 

(top left), 

Lorentz force r Maxc n (top right) and r Total = r advection + r Maxwc n (bottom) in the x-z plane at 
y = and t = 0.5. 





P 


P 


v x 


Vy 


v z 


B x 


By 


B z 


Upstream Region 


0.87 


2.44 








-7.2 


2.87 








Downstream Region Simulation 


1.35 


9.6 





0.29 


-3.58 


5.79 





0.14 


Downstream Region Calculation 


1.81 


12.12 








-3.45 


5.99 









Table III: Values of physical quantities across the MHD shock. All velocities are in the shock rest 
frame, and all magnetic field values are normalized by \/~4m. 



32 



U 



t=0.5 



-5 5 




Figure 19: (color) Contour plot of the specific azimuthal angular momentum prcVip in the x-z plane 
at y = at different times. Left: t = 0. Right: t = 0.5. Net angular momentum is transported 
from the right hand side to the left hand side of the bubble edge. 





t = 


t = 0.25 


t = 0.5 


t = 1.0 


4>t = f BydS 


134.9 


96.9 


88.5 


64.3 



Table IV: Decay of the net toroidal magnetic flux ipt = J B y dS, where only positive B y is selected. 
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